4 Reaction-Diffusion Processes – Advanced Mathematical Biology IV (MATH4411)

Up to this point, we have considered well-mixed chemical reaction processes and this assumption allowed us to neglect any spatial structure in our models. However, many biological phenomena generate spatial pattern and structure; moreover, this spatial structure is often essential to the function of the system. Hence we need to extend our toolkit to include classes of stochastic processes with explicit spatial extent. To this end, we will introduce several stochastic models of the simplest type of spatial movement: diffusion. We then extend our diffusion processes to reaction-diffusion processes in which particles can not only move, but can also react with one another to exhibit more complex dynamics. We illustrate this theory with some applications to pattern formation in biological and ecological systems.

4.1 Modelling diffusion with SDEs

Our first approach to modelling diffusive motion of a particle (or cell, animal, insect, etc.) is to directly model the dynamics of the position of the particle with a stochastic differential equation (SDE). For example, we can denote the position of the particle at time $t$ in 3-dimensional space by $(X(t),Y(t),Z(t))$. If $D>0$ is the diffusion coefficient, the dynamics of the particle are given by

\[ \begin{aligned} X(t+ \Delta t) &= X(t) + \sqrt{2D \Delta t}\, \xi_x, \\ Y(t+ \Delta t) &= Y(t) + \sqrt{2D \Delta t}\, \xi_y, \\ Z(t+ \Delta t) &= Z(t) + \sqrt{2D \Delta t}\, \xi_z, \end{aligned} \tag{4.1}\]

where $\xi_x$, $\xi_y$ and $\xi_z$ are three independent sequences of standard Normal random variables. The $x$, $y$ and $z$ coordinates of the particle are thus fully independent of one another and we are really just considering three separate one dimensional processes; we could construct a diffusion process on $\mathbb{R}^d$ in exactly the same way.

The model of diffusion given by Equation 4.1 is continuous in both space and time (although we choose to represent it using our computational definition). We can thus write down the Fokker-Planck equation(s) for the probability density function for the position of the particle directly from Equation 4.1, i.e.

\[ \begin{aligned} \frac{\partial}{\partial t}P_x(\eta_x,t) &= D \, \frac{\partial^2 }{\partial \eta_x^2}P_x(\eta_x,t), \\ \frac{\partial}{\partial t}P_y(\eta_y,t) &= D \, \frac{\partial^2 }{\partial \eta_y^2}P_y(\eta_y,t), \\ \frac{\partial}{\partial t}P_z(\eta_z,t) &= D \, \frac{\partial^2 }{\partial \eta_z^2}P_z(\eta_z,t), \end{aligned} \tag{4.2}\]

where $P_x(\eta_x,t)$ denotes the marginal density of the $x$ coordinate of the process at time $t$. Since the three components of the process are independent, we can obtain the joint density function of the particle’s position, $P(\eta_x,\eta_y,\eta_z,t)$, by multiplication, i.e.

\[ P(\eta_x,\eta_y,\eta_z,t) = P_x(\eta_x,t)\,P_y(\eta_y,t)\,P_z(\eta_z,t). \]

Hence we can sum the equations in Equation 4.2 to obtain the Fokker-Planck equation for the joint density of the position:

\[ \frac{\partial}{\partial t} P = D\left( \frac{\partial^2}{\partial \eta _x^2} P + \frac{\partial^2}{\partial \eta _y^2} P + \frac{\partial^2}{\partial \eta _z^2} P\right), \]

which is just the heat equation in three-dimensional space. It is straightforward to show that

\[ P(\eta_x,\eta_y,\eta_z,t) = \frac{1}{(4D\pi t)^{3/2}}\exp\left( \frac{-\eta_x^2 - \eta_y ^2 - \eta_z^2 }{4Dt} \right). \tag{4.3}\]

In order to visualize this PDF, we can integrate out the $z$ component to obtain the joint density of the $x$ and $y$ components:

\[ P(\eta_x,\eta_y,t) = \int_{\mathbb{R}}\frac{1}{(4D\pi t)^{3/2}}\exp\left( \frac{-\eta_x^2 - \eta_y ^2 - \eta_z^2 }{4Dt} \right)d\eta_z = \frac{1}{4D\pi t}\exp\left( \frac{-\eta_x^2 - \eta_y ^2}{4Dt} \right). \tag{4.4}\]

Figure 4.1: Comparison of the 2D density Equation 4.4 at time $t = 10$ (left) with the estimated density computed from direct simulations of the process Equation 4.1 (right). Parameters: $D = 0.0001$. [Code: `CH4_SDE_diffusion.m`]

Figure 4.1 shows the two-dimensional density Equation 4.4 at time $t = 10$ (left). From the previous chapter, we know how to simulate the process Equation 4.1 using the computational definition and Figure 4.1 shows the results of estimating the joint density of $(X(t),Y(t))$ via this approach. To estimate the joint density via simulations, we need to discretize space, count the number of sample paths finishing in each box at $t=10$ (see Figure 4.1) and then appropriately normalize these counts for the number of paths and area of the boxes.

We can reduce Equation 4.1 to a one-dimensional model of diffusive motion by just considering the first component, i.e.

\[ X(t+ \Delta t) = X(t) + \sqrt{2D \Delta t}\, \xi. \tag{4.5}\]

The cells or other particles we wish to model often cannot move freely throughout the spatial domain and hence we sometimes need to impose barriers on the domain of $X(t)$ to reflect physical reality. Suppose we want to restrict the domain of $X(t)$ to the interval $[0,L]$ for some $L>0$ by imposing reflecting boundaries at $0$ and $L$. We can simulate sample paths of Equation 4.5 with these two reflecting boundaries using a variation on the SSA introduced in the previous chapter:

At time $t = 0$, set $X(0)$, then:

Step 1: Generate $\xi \sim N(0,1)$.
Step 2: Compute the possible value of $X(t+\Delta t)$:

\[ X(t + \Delta t) = X(t) + \sqrt{2D \, \Delta t}\,\xi, \]
Step 3: If $X(t+\Delta t)\in [0,L]$, then we accept that value.

Otherwise, if $X(t+ \Delta t) < 0$, then set

\[ X(t + \Delta t) = \textcolor{blue}{- X(t)} - \sqrt{2D\,\Delta t}\,\xi. \]

If $X(t+ \Delta t) > L$, then set

\[ X(t + \Delta t) = X(t) + \textcolor{blue}{2(L-X(t))} - \sqrt{2D\,\Delta t}\,\xi. \]

Finally, set $t = t+\Delta t$ and go back to Step 1.

The reflection steps in the SSA above are highlighted in blue with the first instance being a reflection of the process through zero and the second a reflection through the line $\eta_x = L$. Figure 4.2 A shows some sample paths of the process obtained via the SSA above with reflection at $\eta_x = 0$ and $\eta_x = L = 1$.

We can also analyse the 1D diffusion process Equation 4.1 by solving its associated Fokker-Planck equation

\[ \frac{\partial}{\partial t}P_x(\eta_x,t) = D \, \frac{\partial^2 }{\partial \eta_x^2}P_x(\eta_x,t), \quad \text{ for }\eta_x \in (0,L). \]

We additionally need to supply boundary conditions at $\eta_x = 0$ and $\eta_x = L$ to account for the reflection of the process that ensures $X(t) \in [0,L]$ for all $t \geq 0$. From our discussion of reflecting boundary conditions for the Fokker-Planck equation in Chapter 3, we know that we need to specialise condition Equation 3.25; doing so yields the pair of boundary conditions:

\[ \frac{\partial }{\partial \eta_x}P_x(\eta_x,t) \bigg|_{\eta_x = 0} = \frac{\partial }{\partial \eta_x}P_x(\eta_x,t) \bigg|_{\eta_x = L} = 0 \quad \text{for all } t > 0. \tag{4.6}\]

Figure 4.2 A, B and C show how the probability density of $X(t)$, $P(\eta_x,t)$, evolves from an initial distribution that is a Dirac delta at $\eta_x=0.4$. The influence of the reflecting boundaries at $\eta_x=0$ and $\eta_x=1$ become apparent by time $t = 5$ and eventually the distribution of the particle’s position approaches a uniform distribution on $[0,1]$ for large times. The blue histograms in Figure 4.2 are the estimated probability density from the SSA and we observe relatively good agreement between the exact and estimated PDFs at all times ($2,500$ paths). We could improve the quality of the SSA approximation by decreasing the step size in the simulations, at additional computational cost, or by narrowing the width of the histogram bins.

Figure 4.2: A: Sample paths of Equation 4.5 with reflecting boundaries at $0$ and $L = 1$. B/C/D: Solution to the 1D Fokker-Planck equation with reflecting (no-flux) boundary conditions as given by Equation 4.6 at times $t = \{1, \, 5,\, 120\}$. Parameters: $D = 0.0001$. [Codes: `CH4_SDE_reflect.m` and `CH4_FP_reflect.m`]

4.2 A compartmental model of diffusion

In the previous section we modelled diffusive motion by simulating the trajectory of individual particles and then averaging over the trajectories of many particles to understand the properties of the underlying process. We could instead change perspective and consider a collection of particles from the outset and, rather than tracking individual particles, we could record only changes in the number of particles present at each location.

Consider the spatial domain $[0, 1]$ and discretise it into $K = 40$ compartments of equal width, i.e. $h = 1/40$ is our discretisation parameter. Now we can set up a chemical reaction process where each “species” $A_i$ for $i = 1,\dots,40$ counts the number of molecules present in compartment $i$, which is the interval $[(i-1)h, \, ih]$. Next introduce a parameter $d>0$ which we call the “hopping rate” that controls the rate at which particles diffuse into neighbouring compartments. We defer for now the discussion of how to choose the discretisation parameter $h>0$ and the hopping rate $d$ to give a diffusivity rate $D$ corresponding to the diffusion coefficient from the SDE diffusion model Equation 4.1. The dynamics of this process are thus given by:

\[ A_1 \overset{d}{\underset{d}\rightleftarrows} A_2 \overset{d}{\underset{d}\rightleftarrows} \dots\dots \overset{d}{\underset{d}\rightleftarrows} A_{i-1} \overset{d}{\underset{d}\rightleftarrows} A_i \overset{d}{\underset{d}\rightleftarrows} A_{i+1} \overset{d}{\underset{d}\rightleftarrows} \dots\dots \overset{d}{\underset{d}\rightleftarrows} A_{K-1} \overset{d}{\underset{d}\rightleftarrows} A_K. \tag{4.7}\]

This system has $2(K-1)$ reactions but only $K$ distinct propensities. Let $\alpha_i(t) = d A_i(t)$ and note that this propensity represents two reactions:

\[ A_i \xrightarrow{d} A_{i+1} \text{ (jump right)}, \quad \text{and} \quad A_i \xrightarrow{d} A_{i-1} \text{ (jump left)}. \]

Suppose that we start with $N$ molecules in the system, i.e. $\sum_{i=1}^{K} A_i(0) = N$. Then the total propensity of the system is given by

\[ \begin{aligned} \alpha_0(t) &= \sum_{i=1}^{K-1} \alpha_i(t) + \sum_{i=2}^K \alpha_i(t) \\ &= 2 \sum_{i=1}^K \alpha_i(t) - \alpha_1(t) - \alpha_K(t) \\ &= 2d \sum_{i=1}^K A_i(t) - \alpha_1(t) - \alpha_K(t) \\ &= 2 d N - \alpha_1(t) - \alpha_K(t). \end{aligned} \tag{4.8}\]

We can immediately see how to exploit the additional structure of this process to optimize its simulation via the Gillespie algorithm. In particular, the total number of molecules ($N$) is always conserved and only two molecule numbers and two propensity functions change at each reaction. Moreover, as Equation 4.8 illustrates, the total propensity function only changes when there is a reaction impacting compartment $1$ or compartment $K$. Thus we have the following optimized Gillespie algorithm for simulating our compartmental diffusion process:

At time $t = 0$, set the initial molecule numbers in each compartment ($A_i(0)$ for $i=1,\dots,K$), then:

Step 1: Generate $r_1,\,r_2 \sim U([0,1])$.
Step 2: Compute the propensity functions $\alpha_i(t)$ and hence the total propensity function $\alpha_0(t)$ from Equation 4.8.
Step 3: Compute the time of the next reaction, which takes place at time $t+\tau$, by calculating \[ \tau = \frac{1}{\alpha_0(t)}\log(1/r_1). \]
Step 4: If $r_2 < \sum_{i=1}^{K-1} \alpha_i/\alpha_0$, then find $j \in \{1,2,\dots,K-1\}$ such that \[ r_2 \geq \frac{1}{\alpha_0}\sum_{i=1}^{j-1} \alpha_i \text{ and }r_2 < \frac{1}{\alpha_0} \sum_{i=1}^j \alpha_i. \] Then adjust the number of molecules as follows: \[ A_j(t+\tau) = A_j(t)-1, \quad A_{j+1}(t+\tau) = A_{j+1}(t)+1, \] i.e. one molecule in compartment $j$ jumps to the right into compartment $j+1$.
Step 5: If $r_2 \geq \sum_{i=1}^{K-1} \alpha_i/\alpha_0$, then find $j \in \{2,\dots,K\}$ such that \[ r_2 \geq \frac{1}{\alpha_0}\left(\sum_{i=1}^{K-1} \alpha_i + \sum_{i=2}^{j-1} \alpha_i \right) \text{ and }r_2 < \frac{1}{\alpha_0} \left(\sum_{i=1}^{K-1} \alpha_i + \sum_{i=2}^j \alpha_i \right). \] Then adjust the number of molecules as follows: \[ A_j(t+\tau) = A_j(t)-1, \quad A_{j-1}(t+\tau) = A_{j-1}(t)+1, \] i.e. one molecule in compartment $j$ jumps to the left into compartment $j-1$.
Step 6: Set $t = t+\tau$ and go back to Step 1.

Figure 4.3 below shows some simulations of the compartmental diffusion process in panel A. Note that we have not explicitly implemented a reflecting boundary but particles cannot move left from compartment $A_1$ or right from compartment $A_K$ and so the process is naturally confined to the interval $[0,1]$. In panels $B$, $C$ and $D$ of Figure 4.3 we plot the estimated probability density of the process at various times and compare it to the predictions of the corresponding Fokker-Planck equations for the analogous continuum (SDE) model of diffusion. We see excellent agreement between the compartmental and continuum models because we chose the diffusion coefficient $D = d \times h^2$, where $d$ is the compartmental hopping rate and $h$ is our discretisation parameter. But how do we know that this is the correct relationship between the parameters of the discrete and continuous models of diffusion?

To understand how our discrete and continuous models relate to one another, we need to write down the chemical master equations for the compartmental model. To this end, let $p(\textbf{n},t) = \mathbb{P}[\textbf{A}(t) = \textbf{n}]$ denote the joint probability mass function of the compartmental process. Next define the discrete shift operators $R_i,\,L_i: \mathbb{N}^K \mapsto \mathbb{N}^K$ by

\[ \begin{aligned} R_i: [n_1,\dots,n_i,n_{i+1},\dots,n_k] &\mapsto [n_1,\dots,n_i+1,n_{i+1}-1,\dots,n_k], \quad i=1,2,\dots,K-1, \\ L_i: [n_1,\dots,n_{i-1},n_{i},\dots,n_k] &\mapsto [n_1,\dots,n_{i-1}-1,n_{i}+1,\dots,n_k], \quad i=2,3,\dots,K. \end{aligned} \]

The chemical master equations for the process Equation 4.7 are thus given by:

We can then define the mean of the process Equation 4.7 by $\textbf{M}(t) = [M_1(t),\dots,M_K(t)]$ where $M_i(t)$ denotes the mean number of molecules in compartment $i$ and is given by

\[ M_i(t) = \sum_{\textbf{n}} n_i \, p(\textbf{n},t) := \sum_{n_1 = 0}^N \sum_{n_2 = 0}^N \cdots \sum_{n_K = 0}^N n_i \, p(\textbf{n},t). \]

Figure 4.3: A: Sample paths of the compartmental diffusion process with $K = 40$ compartments and $d = D/h^2$. B/C/D: Simulations of the compartmental diffusion process overlaid with solutions to the 1D Fokker-Planck equation with reflecting (no-flux) boundary conditions as given by Equation 4.6 at times $t = \{1, \, 5,\, 120\}$. Parameters: $D = 0.0001$, $h = 0.025$. [Codes: `CH4_compartment_diff.m` and `CH4_compartment_FP.m`]

With some work, we can show that the evolution equations for the mean are given by

\[ \begin{aligned} \frac{d}{dt} M_i(t) &= d\left( M_{i+1} + M_{i-1} - 2M_i \right), \quad i = 2,3,\dots K-1, \\ \frac{d}{dt}M_1 &= d(M_2 - M_1), \quad \frac{d}{dt}M_K = d(M_{K-1} - M_K). \end{aligned} \tag{4.10}\]

To compare this to the continuum model of diffusion studied in the previous section, we need to convert the mean number of molecules, $M_i$, to a concentration. If $c(x,t)$ denotes the concentration of molecules at location $x$, we can make the approximation

\[ c(x_i,t) \approx \frac{M_i(t)}{h} \]

where $x_i$ denotes the center of the $i$th compartment for $i = 1,\dots,K$. Using this approximation, divide Equation 4.10 across by $h$ to obtain

\[ \frac{\partial }{\partial t}c(x_i,t) = d \left( c(x_i + h,t) + c(x_i-h,t) - 2 c(x_i,t) \right). \]

Assuming $c$ is sufficiently smooth, we can Taylor expand the right-hand side to obtain the approximation

\[ \frac{\partial }{\partial t}c(x_i,t) \approx d h^2 \frac{\partial^2 }{\partial x^2} c(x_i,t). \]

Hence if we choose $d = D/h^2$, where $D$ is the continuum diffusion coefficient, we have

\[ \frac{\partial }{\partial t}c(x_i,t) = D \frac{\partial^2 }{\partial x^2} c(x_i,t), \]

and the evolution equation for the molecule concentration will match the corresponding Fokker-Planck equation of the SDE model of diffusion. In other words, the compartmental and continuum diffusion models will agree on average.

To fully characterize the distribution of the compartmental diffusion process we additionally need to understand the variance of the process. The variance vector of the process Equation 4.7 is given by $\textbf{V}(t) = [V_1(t),\dots,V_K(t)]$, where

\[ V_i(t) = \sum_{\textbf{n}}\left(n_i - M_i(t)\right)^2 p(\textbf{n},t) := \sum_{n_1 = 0}^N \sum_{n_2 = 0}^N \cdots \sum_{n_K = 0}^N \left(n_i - M_i(t)\right)^2 \, p(\textbf{n},t). \]

Deriving a system of evolution equations for the variances requires us to define more generally the covariance matrix $V_{i,j}$ by

\[ V_{i,j} = \sum_{\textbf{n}} n_i n_j p(\textbf{n},t) - M_i M_j, \quad i,j = 1,\dots,K. \tag{4.11}\]

The diagonal entries of the covariance matrix Equation 4.11 are the variances defined above. Multiplying the chemical master equations by $n_i^2$ and summing over $\textbf{n}$ yields

\[ \begin{aligned} \frac{d}{dt} \sum_{\textbf{n}} n_i^2 p(\textbf{n},t) &= d \sum_{j = 1}^{K-1} \left( \sum_{\textbf{n}} n_i^2 (n_j+1) p(R_j \textbf{n}) - \sum_{\textbf{n}} n_i^2 n_j p(\textbf{n},t) \right) \\ &\quad+ d \sum_{j = 2}^{K} \left( \sum_{\textbf{n}} n_i^2 (n_j+1) p(L_j \textbf{n}) - \sum_{\textbf{n}} n_i^2 n_j p(\textbf{n},t) \right). \end{aligned} \tag{4.12}\]

Suppose $i \in \{2,\dots,K-1\}$. Evaluate the first term on the right-hand side of Equation 4.12 for $j=i$ and by changing indices $R_i \textbf{n} \mapsto \textbf{n}$:

\[ \begin{aligned} \sum_{\textbf{n}} n_i^2 (n_i+1) p(R_j \textbf{n}) - \sum_{\textbf{n}} n_i^2 n_i p(\textbf{n},t) &= \sum_{\textbf{n}} (n_i-1)^2 (n_i) p(\textbf{n}) - \sum_{\textbf{n}} n_i^2 n_i p(\textbf{n},t) \\ &= \sum_{\textbf{n}} (-2n_i^2 + n_i) p(\textbf{n}) = -2 V_i - 2M_i^2 + M_i. \end{aligned} \]

Now we need to deal with the off-diagonal terms. The term corresponding to $j = i-1$ in the first sum can be rewritten as

\[ \begin{aligned} \sum_{\textbf{n}} n_i^2 (n_{i-1}+1) p(R_{j-1} \textbf{n}) - \sum_{\textbf{n}} n_i^2 n_{i-1} p(\textbf{n},t) &= \sum_{\textbf{n}} (2 n_i n_{i-1} + n_{i-1}) p(\textbf{n}) \\ &= 2 V_{i,i-1} + 2M_i M_{i-1} + M_{i-1}. \end{aligned} \]

Every other term in the first sum on the right-hand side of Equation 4.12 with $j \notin \{i,i-1\}$ is zero. The second sum on the right-hand side of Equation 4.12 is handled analogously by addressing first the cases $j = i$ and $j = i+1$, and showing that all other terms are zero. Carrying out the requisite algebra yields

\[ \begin{aligned} \frac{d}{dt} \sum_{\textbf{n}} n_i^2 p(\textbf{n},t) &= d \left( 2V_{i,i-1} + 2M_i M_{i-1} + M_{i-1} - 2V_i - 2M_i^2 + M_i \right) \\ &\quad + d\left( 2V_{i,i+1} + 2M_i M_{i+1} + M_{i+1} - 2V_i - 2M_i^2 + M_i \right). \end{aligned} \]

Finally, we arrive at the following set of equations for the variances of the compartments:

\[ \begin{aligned} \frac{d}{dt} V_i &= 2d\left( V_{i,i+1} + V_{i,i-1} - 2V_i \right) + d\left( M_{i+1} + M_{i-1} + 2M_i \right), \quad i = 2,\dots,K-1,\\ \frac{d}{dt}V_1 &= 2d\left( V_{1,2} - V_1 \right) + d\left( M_{2} + M_1 \right), \\ \frac{d}{dt}V_K &= 2d\left( V_{K,K-1} - V_K \right) + d\left( M_{K-1} + M_K \right), \end{aligned} \tag{4.13}\]

However, the system Equation 4.13 is not a fully closed set of equations because it involves covariance terms of the form $V_{i,j}$ for $i \neq j$. Hence we need to also write down evolution equations for the off-diagonal covariance matrix terms; this can be done by repeating variations on the arguments above.

We can also consider the steady state versions of the mean and variance equations to study the long-term behaviour of the compartmental diffusion process. Doing so yields

\[ \bar{M}_i = \frac{N}{K} \quad \text{for } i = 1,\dots,K, \] and \[ \bar{V}_i = \frac{N}{K} - \frac{N}{K^2}, \quad \bar{V}_{i,j} = \frac{-N}{K^2}, \quad i \neq j. \]

In fact, solving the chemical master equations directly shows that the stationary distribution of the compartmental diffusion process is a multinomial distribution, i.e. \[ p_s(\textbf{n}) = \frac{C}{n_1 ! n_2 ! \cdots n_K !}, \] for an appropriate normalising constant $C > 0$.

4.3 SSAs for reaction-diffusion processes

4.3.1 Compartmental models

We begin with a discrete-space approach that extends our compartmental model of diffusion to incorporate simple reaction dynamics (zero and first order reactions). We begin by defining an essentially one-dimensional spatial domain

\[ \Omega = [0,L] \times [0,h] \times [0,h] \]

where we assume $h \ll L$ so that spatial dynamics are only occurring in the $x$ dimension. Formally, we can assume that we have periodic boundary conditions in the $y$ and $z$ directions. Next divide the interval $[0,L]$ into $K$ compartments of width $h$, as we did before, and let $A_i(t)$ denote the number of molecules present in compartment $i$ at time $t$. Note that compartment $i$ is the box $[(i-1)h, ih] \times [0,h] \times [0,h]$ and hence the volume of compartment $i$ is $h^3$. The dynamics of the process will be as follows:

\[ \begin{aligned} A_1 \overset{d}{\underset{d}\rightleftarrows} & A_2 \overset{d}{\underset{d}\rightleftarrows} \dots\dots \overset{d}{\underset{d}\rightleftarrows} A_{i-1} \overset{d}{\underset{d}\rightleftarrows} A_i \overset{d}{\underset{d}\rightleftarrows} A_{i+1} \overset{d}{\underset{d}\rightleftarrows} \dots\dots \overset{d}{\underset{d}\rightleftarrows} A_{K-1} \overset{d}{\underset{d}\rightleftarrows} A_K, \\ & A_i \xrightarrow{k_1} \emptyset, \qquad i = 1,\dots,K, \quad \emptyset \xrightarrow{k_2} A_i, \quad i = 1, \dots,K/5. \end{aligned} \tag{4.14}\]

To summarize the dynamics above:

Molecules diffuse at rate $d$ between all compartments (and we can choose $d = D/h^2$ in order to relate this model to a continuum diffusion approach),
molecules are produced in the subdomain $[0,L/5]$ at rate $k_2$ per compartment,
all molecules in all compartments are subject to degradation at rate $k_1$.

In order to analyse the dynamics of the process defined by the system of reactions Equation 4.14, we need to write down the reaction-diffusion (RD) master equations for the process. The RD master equations are:

\[ \begin{aligned} \frac{d}{dt}p(\textbf{n}) &= d \sum_{j=1}^{K-1} \left( (n_j + 1)p(R_j \textbf{n}) - n_j p(\textbf{n}) \right) \\ &\quad + d \sum_{j=2}^{K} \left( (n_j + 1)p(L_j \textbf{n}) - n_j p(\textbf{n}) \right) \\ &\quad +k_1 \sum_{i = 1}^K \left( (n_i + 1)p(n_1,\dots,n_i + 1,\dots,n_k) - n_i p(\textbf{n}) \right) \\ &\quad + k_2 h^3 \sum_{i = 1}^{K/5} \left( p(n_1,\dots,n_i - 1,\dots,n_k) p(\textbf{n}) \right). \end{aligned} \tag{4.15}\]

The first two terms on the right-hand side of Equation 4.15 represent diffusive reactions, the third term represents the degradation events and the final term captures the impact of production reactions in the subdomain $[0,L/5]$. We can once more define the mean vector $\textbf{M}(t) = [M_1(t),\dots,M_K(t)]$ and write down a set of evolution equations for the components of $\textbf{M}(t)$. Using the fact that all reactions in Equation 4.14 are either zero or first order, we can directly write down the evolution equations for the stochastic means:

\[ \begin{aligned} \frac{d}{dt}M_1 &= d(M_2 - M_1) + k_2 h^3 - k_1 M_1,\\ \frac{d}{dt} M_i &= d\left( M_{i+1} + M_{i-1} - 2M_i \right) + k_2 h^3 - k_1 M_i, \quad i = 2,3,\dots K/5, \\ \frac{d}{dt} M_i &= d\left( M_{i+1} + M_{i-1} - 2M_i \right) - k_1 M_i, \quad i = K/5 +1,\dots K-1, \\ \frac{d}{dt}M_K &= d(M_{K-1} - M_K) - k_1 M_K. \end{aligned} \]

To better compare the stochastic means with a spatially continuous approximation, we let $\tilde{M}_i(t) = M_i(t)/h^3$ to move from molecule numbers to molecule concentrations. This volumetric scaling and choosing $d = D/h^2$ yields:

\[ \begin{aligned} \frac{d}{dt}\tilde{M}_1 &= D \frac{\tilde{M}_2 - \tilde{M}_1}{h^2} + k_2 - k_1 \tilde{M}_1,\\ \frac{d}{dt} \tilde{M}_i &= D \frac{\tilde{M}_{i+1} + \tilde{M}_{i-1} - 2\tilde{M}_i }{h^2} + k_2 - k_1 \tilde{M}_i, \quad i = 2,3,\dots K/5, \\ \frac{d}{dt} \tilde{M}_i &= D \frac{\tilde{M}_{i+1} + \tilde{M}_{i-1} - 2\tilde{M}_i}{h^2} - k_1 \tilde{M}_i, \quad i = K/5 +1,\dots K-1, \\ \frac{d}{dt} \tilde{M}_K &= D \frac{\tilde{M}_{K-1}-\tilde{M}_K}{h^2} - k_1 \tilde{M}_K. \end{aligned} \tag{4.16}\]

From the definition of the derivative, we have the following first-order finite-difference approximation:

\[ \frac{d}{dx}f(x) = \lim_{h \downarrow 0}\frac{f(x+h)-f(x)}{h} \approx \Delta f(x) := \frac{f(x+h) - f(x)}{h}. \]

Hence we can make the same finite-difference approximation of the second derivative:

\[ \frac{d^2}{dx^2}f(x) \approx \Delta^2 f(x) = \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}. \]

Therefore we can interpret the evolution equations Equation 4.16 as a finite-difference approximation of the following (deterministic) reaction-diffusion PDE:

\[ \frac{\partial }{\partial t}m(x,t) = D\, \frac{\partial^2}{\partial x^2}m(x,t) + k_2 \mathbf{1}_{\{x \in [0,L/5]\}} - k_1 m(x,t), \quad x \in (0,L), \tag{4.17}\]

with reflecting boundary conditions:

\[ \frac{\partial}{\partial x}m(x,t) {\bigg \vert}_{x = 0} = \frac{\partial}{\partial x}m(x,t) {\bigg \vert}_{x = L} = 0, \] where $\mathbf{1}$ denotes the indicator function.

Figure 4.4 below shows the results of a single simulation of the reaction-diffusion process Equation 4.14 with the estimated PMF plotted at different times in blue. The solution to the approximating reaction-diffusion PDE Equation 4.17 is plotted for comparison in red and we observe quite close agreement between the two models, even for a single realisation of the stochastic process. Of course, the reaction-diffusion approximation is only expected to capture the mean behaviour and will not model fluctuations about the mean of the process.

Figure 4.4: Simulations of the compartmental reaction-diffusion process Equation 4.14 with the estimated PMF in blue and the PDE approximation Equation 4.17 shown in red. Parameters: $D = 0.0001$, $h = 0.01$, $k_1 = 0.001$, $k_2 = 0.3125$, $L=1$. [Code: `CH4_compartment_RD.m`]

4.3.2 SDE-based reaction-diffusion models

As with the modelling of diffusion, we could take the alternative approach of tracking the movements and reactions of individual particles, rather than working with a discretised spatial domain of compartments. To illustrate this alternative approach, we will build a spatially continuous model of the compartmental reaction-diffusion process Equation 4.14. We will once more consider the spatial domain to be $\Omega = [0,L] \times [0,h] \times [0,h]$ with $h \ll L$ so that we only model spatial movement in the $x$ dimension. We now want the continuum model to capture the following attributes:

Molecules diffuse at rate $D$ on $[0,L]$ with reflecting boundaries at $x = 0$ and $x = L$,
molecules are produced in the subdomain $[0,L/5]$ at rate $k_2$,
all molecules are subject to degradation at rate $k_1$.

To simulate the stochastic process outlined above, we will discretise time (rather than space) by introducing a finite-size time step $\Delta t>0$. The spatial dynamics of each particle $X_i$ will be given by the following diffusive SDE:

\[ X_i(t + \Delta t) = X_i(t) + \sqrt{2 D \Delta t}\, \xi. \tag{4.18}\]

Then at each time $t$ we will evolve the process forward in time by carrying out the following steps:

Step 1: For each molecule, compute its $x$ position at time $t+\Delta t$ according to Equation 4.18 (obeying the reflecting boundary conditions at $x = 0$ and $x = L$).
Step 2: For each molecule, generate $r_1 \sim U([0,1])$. If $r_1 < k_1 \Delta t$, remove the molecule (degradation reaction occurs).
Step 3: Generate $r_2 \sim U([0,1])$. If $r_2 < (k_2 h^2 L/5)\Delta t$, then generate $r_3 \sim U([0,1])$ and introduce a new molecule at position $x = r_3 L / 5$ (production reaction occurs).

We don’t further specify the details of step 1, since we have previously detailed the algorithm to simulate the diffusive SDE Equation 4.18.

If we simulate the reaction-diffusion process using this algorithm and estimate the probability density of the number of molecules we obtain results virtually indistinguishable to those shown in Figure 4.4 (assuming we choose $D = dh^2$). But why should this new algorithm produce statistically accurate realisations of the diffusion-production-degradation process?

The algorithm above essentially couples our computational definition of SDEs to the “naive” SSA for chemical reaction systems that we introduced at the very beginning of the course. If we choose $\Delta t$ sufficiently small that $\Delta t k_1 \ll 1$, then we correctly reproduce the degradation reaction in step 2, because

\[ \mathbb{P}\left[\text{molecule $i$ gets degraded in }[t,\,t+\Delta t] \right] \approx k_1 \Delta t + O(\Delta t^2). \]

Similarly, new molecules are produced in the region $[0,L/5]$, which has volume $h^2 L/5$, and so

\[ \mathbb{P}\left[\text{new molecule produced in }[t,\,t+\Delta t] \right] \approx k_2 h^2 L/5 \Delta t + O(\Delta t^2). \]

As with the degradation reactions, we will need to choose the time step sufficiently small to ensure accuracy of the algorithm. In this case, we need $k_2 h^2 L/5 \Delta t \ll 1$ and this can also be achieved by choosing $h$ sufficiently small since $h$ is the spatial discretisation parameter. Finally, if a molecule is produced at a given time step, then it was equally likely to have been produced at any location in $[0,L/5]$ and so we simply assign it to a uniformly random location in $[0,L/5]$ in the second part of step 3.

Note that the ordering of the steps in the algorithm here does not matter. If we move molecule $i$ in the SDE updating step, but then that molecule gets degraded at the same time step, it does not impact the dynamics of any other molecule. We also check the degradation reaction for every molecule at each time step so it is possible for multiple molecules to be degraded at each time step but this is not a problem since the molecules are not interacting with each other! In other words, we don’t have to worry about a case in which two molecules get degraded in one time step but the presence of one of these molecules would have influenced the subsequent dynamics of the other; we will come back to these potential issues in the next section.

One other difference to this continuum space approach in which we track molecule trajectories is that we need to dynamically track the number of particles. From a computational perspective, the size of our system is not fixed and we need to dynamically update the number of SDEs that we are solving and discontinue any SDEs that correspond to molecules that have been degraded. We also need to introduce and solve a new SDE when a new molecule enters the system via a production reaction. There are numerous ways to handle this complication during the implementation of the algorithm but it is worth factoring in when comparing the computational cost of alternative modelling approaches.

4.3.3 Reaction-Diffusion processes with higher-order reactions

In the preceding sections, we introduced a discrete space and a continuous space approach to analysing and simulating stochastic reaction-diffusion processes. In both cases, we restricted our examples to zero and first order reactions, and both approaches were straightforward to implement. Of course, most models of interest in biology will involve higher-order reactions (i.e. interactions between at least two different species) and this forces us to answer some questions which simply did not arise with only zero or first order reactions. In particular, if two molecules are required for a reaction to occur, how close do they need to be in space for the reaction to occur? Should the reaction probability depend on the distance between molecules or always occur if they are “close enough”?

Consider the following illustrative example of a compartmental reaction-diffusion process: There are two types of molecules, type A and type B, present in the pseudo-one-dimensional domain

\[ \Omega = [0,L] \times [0,h] \times [0,h], \]

and as usual we only model their movement in the $x$ dimension with the assumption that $h \ll L$. The domain $[0,L]$ is discretised into $K$ compartments of width $h$. The number of $A$ molecules in compartment $i$, the interval $[(i-1)h,ih]$, at time $t$ is given by “species” $A_i(t)$ and, similarly, $B_i(t)$ denotes the number of $B$ molecules in compartment $i$. In terms of reactions, posit that the process has the following dynamics:

$A$ molecules move to adjacent compartment at hopping rate $d_A$ ($d_A = D_A/h^2$) and $B$ molecules have hopping rate $d_B$ ($d_B = D_B/h^2$),
molecules can only react with other molecules in their own compartment,
Two $A$ molecules react at rate $k_1$ to produce a product not of interest,
An $A$ and a $B$ molecule react at rate $k_2$ to produce a product not of interest,
$A$ molecules are degraded at rate $k_3$ in $[0,L]$,
$B$ molecules are degraded at rate $k_4$ in $[0,L]$,
$A$ molecules are produced at rate $k_5$ in $[0,L]$,
$B$ molecules are produced at rate $k_6$ in $[3L/5,L]$,

The second bullet point above is making explicit the crucial assumption about which molecules can react with each other, i.e. only those within the same spatial compartment. This is essentially assuming that molecules are well-mixed within each compartment, so that the process has two different spatial scales. The macro scale is the interval $[0,L]$ and the molecules are not well-mixed on this scale (we can have spatial pattern), but on the micro scale (within a compartment) molecules are well-mixed spatially and there is no spatial structure within the compartments.

In our usual reaction notation, we can represent this process as follows:

\[ \begin{aligned} A_1 \overset{d_A}{\underset{d_A}\rightleftarrows} & A_2 \overset{d_A}{\underset{d_A}\rightleftarrows} \dots\dots \overset{d_A}{\underset{d_A}\rightleftarrows} A_{i-1} \overset{d_A}{\underset{d_A}\rightleftarrows} A_i \overset{d_A}{\underset{d_A}\rightleftarrows} A_{i+1} \overset{d_A}{\underset{d_A}\rightleftarrows} \dots\dots \overset{d_A}{\underset{d_A}\rightleftarrows} A_{K-1} \overset{d_A}{\underset{d_A}\rightleftarrows} A_K, \\ B_1 \overset{d_B}{\underset{d_B}\rightleftarrows} & B_2 \overset{d_B}{\underset{d_B}\rightleftarrows} \dots\dots \overset{d_B}{\underset{d_B}\rightleftarrows} B_{i-1} \overset{d_B}{\underset{d_B}\rightleftarrows} B_i \overset{d_B}{\underset{d_B}\rightleftarrows} B_{i+1} \overset{d_B}{\underset{d_B}\rightleftarrows} \dots\dots \overset{d_B}{\underset{d_B}\rightleftarrows} B_{K-1} \overset{d_B}{\underset{d_B}\rightleftarrows} B_K, \\ & A_i + A_i \xrightarrow{k_1} \emptyset, \quad A_i + B_i \xrightarrow{k_2} \emptyset, \quad i = 1,2,\dots,K, \\ & A_i \xrightarrow{k_3} \emptyset, \quad B_i \xrightarrow{k_4} \emptyset, \quad \emptyset \xrightarrow{k_5} A_i, \quad i = 1,\dots,K, \\ & \emptyset \xrightarrow{k_6} B_i, \quad i = 3K/5 + 1, \dots,K. \end{aligned} \tag{4.19}\]

We now have a (potentially very high-dimensional) reaction-diffusion process that can be simulated using the standard Gillespie algorithm. Since we have second-order reactions, we can no longer easily write down a closed system of equations for the stochastic mean of the process. We can however make a law of mass action approximation of the process based on an extension of the mass action principle we used for well-mixed systems. Firstly, let $a(x,t)$ and $b(x,t)$ denote the concentrations of molecules of species $A$ and species $B$ at location $x \in [0,L]$. These concentrations can be approximated from the compartmental model as

\[ a(x_i,t) \approx \frac{A_i(t)}{h^3}, \quad b(x_i,t) \approx \frac{B_i(t)}{h^3}, \]

where $x_i$ denotes the center of compartment $i$. If we set $d_A = d_B = 0$, then the standard law of mass action would give us the following approximate model for the average number of particles at each location:

\[ \begin{aligned} \frac{\partial}{\partial t}a(x,t) &= -2 k_1 a^2 - k_2 a b - k_3 a + k_5, \\ \frac{\partial}{\partial t}b(x,t) &= -k_2 a b + k_6\mathbf{1}_{\{x \in [3L/5,L]\}} - k_4 b, \end{aligned} \tag{4.20}\]

for each $x\in [0,L]$. To extend the law of mass action to our spatial model, we simply include diffusion of both species $A$ and $B$ at the appropriate rates, and add the no flux boundary conditions at $x = 0$ and $x = L$, to yield:

\[ \begin{aligned} \frac{\partial}{\partial t}a(x,t) &= -2 k_1 a^2 - k_2 a b - k_3 a + k_5 + D_A \frac{\partial^2}{\partial x^2} a(x,t), \\ \frac{\partial}{\partial t}b(x,t) &= -k_2 a b + k_6\mathbf{1}_{\{x \in [3L/5,L]\}} - k_4 b + D_B \frac{\partial^2}{\partial x^2} b(x,t),\\ \frac{\partial }{\partial x}a(x,t) {\bigg \vert}_{x = 0} &= \frac{\partial }{\partial x}a(x,t) {\bigg \vert}_{x = L} = \frac{\partial }{\partial x}b(x,t) {\bigg \vert}_{x = 0} = \frac{\partial }{\partial x}b(x,t) {\bigg \vert}_{x = L} = 0. \end{aligned} \tag{4.21}\]

Figure 4.5 below shows the results of simulating one realisation of the stochastic process and estimating the PMF at time $t = 200$, along with the solution of the associated law of mass action approximation from Equation 4.21. Because we are comparing a model of the approximate mean behaviour with a single realisation of the process, there is some disagreement due to the noisy nature of the underlying process but both solutions qualitatively agree in terms of the respective domains of dominance of the $A$ and $B$ molecules.

Figure 4.5: Simulations of the compartmental reaction-diffusion process Equation 4.19 with the estimated PMF in blue and the PDE approximation Equation 4.17 shown in red. Parameters: $D_A = D_B = 0.0001$, $h = 0.01$, $k_1 = 0.03$, $k_2 = 0.3$, $k_3 = 0.0001$, $k_4 = 0.0001$, $k_5 = 0.0000001$, $k_6 = 0.000001$, $L=1$. [Code: `CH4_RD_compartment.m`]

To conclude our discussion of the compartmental approach to modelling reaction-diffusion processes, we need to address a key limitation of the method outlined above. Specifically, we need to discuss the accuracy of the method as it relates to the choice of the compartment size $h$ (also thought of as the spatial discretisation parameter). We saw previously that we could derive reaction-diffusion PDEs in the limit as $h \downarrow 0$ when we have only zero and first order reactions. In fact, in these cases, the reaction-diffusion PDEs are exactly the evolution equations for the stochastic mean of the system as $h \downarrow 0$. However, when we have higher-order (nonlinear) reactions, there is trade-off between accurately modelling the diffusion of molecules and the reactions of molecules as $h \downarrow 0$. In particular, we model diffusion ever more accurately as we decrease $h$, but we model nonlinear reactions less accurately as we decrease $h$.

To understand this drawback of compartmental diffusion, we will revisit the stochastic dimerisation process from Chapter 1, i.e. the well-mixed single-species chemical reaction process with dynamics:

\[ A + A \xrightarrow{k_1} \emptyset, \quad \emptyset \xrightarrow{k_2} A. \tag{4.22}\]

In Chapter 1 we saw that we could express the stationary distribution of the process Equation 4.22 in terms of the modified Bessel function of the first kind and we hence calculated the stationary mean of the process, $M_s$, exactly. This process can thus serve as a good benchmark since we understand its properties fully. We now extend the dimerisation process to the spatial domain

\[ \Omega = [0,L] \times [0,L] \times [0,L], \]

which we discretise into $K$ compartments of width $h$ in each dimension. In other words, we divide $\Omega$ into $K^3$ cubes each of volume $h^3$, where compartment $(i,j,k)$ is the compact interval $[(i-1)h,ih]\times[(j-1)h,jh]\times[(k-1)h,kh]$ for $(i,j,k) \in \{1,\dots,K\}^3$. We then let $A_{i,j,k}(t)$ denote the number of molecules of species $A$ present in the compartment $(i,j,k)$ at time $t$. As in the previous example, only molecules in the same compartment can react with one another so the reactions in Equation 4.22 become:

\[ A_{i,j,k} + A_{i,j,k} \xrightarrow{k_1} \emptyset, \quad \emptyset \xrightarrow{k_2} A_{i,j,k}. \tag{4.23}\]

We also assume that molecules can diffuse at rate $D_A$ ($= d_A h^2$), but since there is production and degradation across the entire domain we don’t expect this to result in any spatial structure in the solution, other than that resulting from stochastic fluctuations. In three dimensions, diffusion corresponds to 6 possible movements of a molecule; the cube of volume $h^3$ that the molecule resides in has 6 faces (unless on a boundary) and a molecule may diffuse into a neighbouring cube corresponding to moving through each one of these faces.

The propensity functions for degradation and production in each compartment of our spatial process are:

\[ \alpha_{i,j,k,1}(t) = \frac{k_1 A_{i,j,k}(t) (A_{i,j,k}(t)-1)}{h^3}, \quad \alpha_{i,j,k,2}(t) = k_2 h^3, \]

and the propensity function for diffusion in compartment $(i,j,k)$ is $D_A A_{i,j,k}(t)/h^2 = d_A A_{i,j,k}(t)$. We can also compute the total number of $A$ molecules in the system at time $t$ as

\[ A(t) = \sum_{i,j,k} A_{i,j,k}(t) \]

and the long-run stochastic mean can then be defined as

\[ M(h) := \lim_{t \to \infty} \mathbb{E}[A(t)], \]

where we write the long-run mean as a function of $h$ in anticipation of studying its behaviour as $h \downarrow 0$. If our SSA works as we would hope, we should obtain $\lim_{h \to 0}M(h) = M_s$.

Figure 4.6: Simulations of the compartmental spatial dimerisation process using the standard Gillespie algorithm ($\beta = 0$) and the modified algorithm with rescaled dimerisation propensities ($\beta = 0.275$). Parameters: $D_A = 0.0001$, $L=1$.

We can simulate the spatial dimerisation process via the Gillespie algorithm by choosing $h = L/K$ and increasing the value of $K$ to send $h \to 0$. The results of these simulations are plotted in Figure 4.6 with the red dots indicating the estimated value of $M(h)$ from standard Gillespie simulations and the solid magenta line showing the true value of $M_s$ for the well-mixed version of the dimerisation process. Since we have only a single species of molecules diffusing, we do not expect any spatial structure or influence on the total molecule count and $M_s$ and $M(h)$ should agree, but we actually see $M(h)$ increasing and diverging from $M_s$ as $h$ decreases. It turns out that $\lim_{h \to 0}M(h)$ does not exist and the SSA is not convergent as $h \to 0$. This is because the dimerisation reaction doesn’t occur often enough as $h \to 0$ and hence we conclude that the compartmental SSA approach is only valid for $h$ sufficiently large! This is exactly the opposite situation when we use a finite difference scheme to numerically approximate the solution to a PDE; we expect to see convergence of the scheme as we refine the mesh by taking the spatial discretisation parameter $h$ smaller and smaller.

The fact that we cannot take $h$ arbitrarily small is immediately disconcerting since we want $h$ small to ensure that we model diffusion accurately. We generally want $h \ll L$ to be able to accurately reflect spatial variation in the solution. However, for this dimerisation process we need to take $h \gg k_1 / D_A$ in order to ensure that the nonlinear reactions aren’t artificially suppressed, leading to a range of acceptable values for $h$:

\[ k_1 / D_A \ll h \ll L. \]

Similarly, if we instead had a two species reaction of the form:

\[ A + B \xrightarrow{k} C, \]

with $A$ and $B$ diffusing at rates $D_A$ and $D_B$ respectively, then we would need $h \gg k/(D_A+D_B)$ to accurately reflect this reaction. Erban and Chapman (2009) have shown that for the spatial dimerisation process, it is possible to remedy this issue by rescaling the propensity function of the dimerisation reaction to

\[ \alpha_{i,j,k,1}(t) = A_{i,j,k}(t) (A_{i,j,k}(t) - 1)\frac{D_A k_1 }{D_A h^3 - \beta k_1 h^2}. \]

For $\beta = 0$, this is the standard Gillespie algorithm but for $\beta > 0$ the algorithm can be refined to accurately recover the value of $M_s$, as shown by the black dots in Figure 4.6. For a given value of $K$, the optimal value of $\beta$ can be determined analytically for this process. However, this is a challenging problem in general, and the optimal choice of $h$ for accuracy in both the reaction and diffusive dynamics is an area of active research.

4.4 Applications to biological pattern formation

There are many examples of regular spatial pattern formation in biology and ecology, including zebra stripes, spotted and striped fish, and vegetation patterns. Turing’s pattern forming mechanism has long been studied as a simple and parsimonious way for self-organized patterns to form in the absence of external cues or other exogenous spatial structure. Turing patterns are typically studied mathematically in PDE-based models, but, as we have seen earlier in this chapter, these models are often coarse-grained approximations that capture average behaviour and neglect fluctuations and finite-size effects. To conclude our investigations of spatially extended stochastic processes, we will thus present an example of a “Turing pattern” in a stochastic reaction-diffusion system.

Consider two species of molecules, $A$ and $B$, reacting and diffusing in our standard pseudo one-dimensional domain

\[ \Omega = [0,L] \times [0,h] \times [0,h], \]

with the assumption that $h \ll L$. We only model diffusive movement in the $x$ coordinate and discretise the interval $[0,L]$ into $K$ compartments, with compartment $i$ given by $[(i-1)h,ih]$. We choose Schnakenberg reaction kinetics so that the dynamics of the model are given by:

\[ \begin{aligned} A_1 \overset{d_A}{\underset{d_A}\rightleftarrows} & A_2 \overset{d_A}{\underset{d_A}\rightleftarrows} \dots\dots \overset{d_A}{\underset{d_A}\rightleftarrows} A_{i-1} \overset{d_A}{\underset{d_A}\rightleftarrows} A_i \overset{d_A}{\underset{d_A}\rightleftarrows} A_{i+1} \overset{d_A}{\underset{d_A}\rightleftarrows} \dots\dots \overset{d_A}{\underset{d_A}\rightleftarrows} A_{K-1} \overset{d_A}{\underset{d_A}\rightleftarrows} A_K, \\ B_1 \overset{d_B}{\underset{d_B}\rightleftarrows} & B_2 \overset{d_B}{\underset{d_B}\rightleftarrows} \dots\dots \overset{d_B}{\underset{d_B}\rightleftarrows} B_{i-1} \overset{d_B}{\underset{d_B}\rightleftarrows} B_i \overset{d_B}{\underset{d_B}\rightleftarrows} B_{i+1} \overset{d_B}{\underset{d_B}\rightleftarrows} \dots\dots \overset{d_B}{\underset{d_B}\rightleftarrows} B_{K-1} \overset{d_B}{\underset{d_B}\rightleftarrows} B_K, \\ & 2A_i + B_i \xrightarrow{k_1} 3A_i, \quad \emptyset \xrightarrow{k_2} A_i, \quad i = 1,2,\dots,K, \\ & A_i \xrightarrow{k_3} \emptyset, \quad \emptyset \xrightarrow{k_4} B_i, \quad i = 1,\dots,K. \end{aligned} \tag{4.24}\]

When choosing the hopping rates, $d_A$ and $d_B$, we employ the standard relationship between the compartmental hopping rates and the continuum diffusion coefficients, i.e.

\[ d_A = \frac{D_A}{h^2}, \quad d_B = \frac{D_B}{h^2}. \]

In the standard nomenclature of Turing pattern formation, species $A$ should act as the slower diffusing activator and species $B$ should play the role of the more quickly diffusing inhibitor species. Thus we choose the ratio of the diffusion coefficients as

\[ \frac{D_B}{D_A} = 100, \quad \text{with }D_A = 10^{-5} \text{ and }D_B = 10^{-3}. \tag{4.25}\]

For this example, the reaction rates are chosen according to:

\[ k_1/h^6 = 10^{-6}, \quad k_2 h^3 = 1, \quad k_3 = 0.02, \quad k_4 h^3 = 3. \]

Figure 4.7 below shows a single realisation of the process Equation 4.24 with $K = 40$ compartments and the reaction and diffusion rates given above ($L = 1$, $h = 0.025$). Panels A and B show the time-space evolution of the process from a homogeneous initial condition. After around 500 seconds, we begin to see the emergence of spatial structure in the abundances of both the $A$ and $B$ molecules, although there is significant variation over time due to the underlying stochasticity of the process. Panels C and D show the state of the system at time $t = 2000$ seconds and we observe clear spatial patterning, especially in species $A$, which has very large amplitude patterns. Since this is a single realisation of the process, the roughness of the patterns is to be expected and we could average over many realisations if we wanted to estimate the mean behaviour, which would display much greater regularity.

Figure 4.7: A/B: Space-time plots of the dynamics of the process Equation 4.24 (single realisation). C/D: Plots of the distribution of $A$ and $B$ molecules at time $t=2000$ (in seconds). [Code: `CH4_Schnakenberg.m`]

To gain additional insight into the process Equation 4.24, we can write down a relatively simple mass action approximation of the average behaviour in the usual way. Letting $a(x,t)$ and $b(x,t)$ denote the $A$ and $B$ molecules concentrations, our mass action model is:

\[ \begin{aligned} \frac{\partial}{\partial t} a &= k_1 a^2 b + k_2 - k_3 a + D_A \frac{\partial^2 }{\partial x^2} a, \quad x \in (0,L),\\ \frac{\partial}{\partial t} b &= -k_1 a^2 b + k_4 + D_B \frac{\partial^2 }{\partial x^2} b, \qquad x \in (0,L),\\ \frac{\partial}{\partial x} a(x,t) {\big \vert}_{x= 0} &= \frac{\partial}{\partial x} a(x,t) {\big \vert}_{x=L} = \frac{\partial}{\partial x} b(x,t) {\big \vert}_{x= 0} = \frac{\partial}{\partial x} b(x,t) {\big \vert}_{x= L} = 0. \end{aligned} \tag{4.26}\]

Linear stability analysis of the approximate deterministic model Equation 4.26 shows that, for our chosen parameter set, it has a spatially homogeneous equilibrium solution corresponding to

\[ \bar{A} = 200, \quad \bar{B} = 75, \]

molecules per volume $h^3$. The spatially homogeneous solution is stable for $D_A = D_B = 0$ but unstable for the diffusion coefficient values in Equation 4.25. Hence the results shown in Figure 4.7 really do represent a classical example of a spatially homogeneous equilibrium that is destabilised by a diffusion driven instability.

Exercise 4.1

Solve the mass action PDE Equation 4.26 numerically using this VisualPDE link. Do the PDE dynamics qualitatively agree with those of the stochastic diffusion model?
Open the course Github page and use the MATLAB script CH4_Schnakenberg.m to compare the dynamics of the stochastic compartmental diffusion process with the dynamics of the PDE Equation 4.26 in VisualPDE for a range of parameter values.
Can you numerically identify the smallest value of $D_A$ at which we no longer observe pattern formation in each model?

Bonus: Can you analytically calculate the value of $D_A$ at which the homogeneous solution of the PDE Equation 4.26 becomes unstable?

Note: In VisualPDE notation, the species $A$ is modelled as $u$ and $B$ is modelled as $v$. Similarly, the diffusion coefficients $D_A$ and $D_B$ are written as $D_u$ and $D_v$ in VisualPDE.

Knowledge checklist

Key topics:

Diffusion based on SDEs (SSAs and exact solutions)
Compartmental diffusion (SSAs and CME analysis)
Reaction-Diffusion processes (SSAs, higher-order reactions, tradeoffs)
Mass action approximations of reaction-diffusion processes.

Key skills:

Represent spatial movement processes as chemical reaction processes.
Use the Chemical Master Equations to analyse reaction-diffusion processes:
- computing stationary distributions from CMEs.
- calculating and solving moment equations.
Write and analyse SSAs, and propose alternative SSAs for reaction-diffusion processes (e.g. Q3 on Problem Sheet 4).
Write and analyse the mass-action (PDE) approximation for a given reaction-diffusion process.
Discuss advantages and disadvantages of different approaches to modelling diffusion as a stochastic process (considering modelling implications and computational costs).